/*
 * main.cpp
 *
 *  Created on: Jan 7, 2011
 *      Author: gibson
 */

#include "InlineStream.h"
#include "PeakReference.h"
#include "noiseReference.h"
#include "averageTwo.h"
#include <iostream>
#include <fstream>
#include <cmath>
#include "roughGaussian.h"
#include "TF1.h"
#include "TH1I.h"
#include "TFitResultPtr.h"
#include "TFitResult.h"
#include "TApplication.h"
#include "TCanvas.h"
#include <time.h>
#include "TGraph.h"
#include "TFile.h"
#include "TMultiGraph.h"

using namespace std;

void addToList(double**& list, int& count, double time, double height,
		double chi2);
void RemoveFromList(double**& list, int& countOfList, int indexToRemove);
void EvaluateAccuracy(double** exactPeakInfo, int exactPeakCount,
		double**& fittedPeakInfo, int& FittedPeakCount);
void RootFineFit(TH1I* histo, int numBins, int roughPosition,
		double roughHeight, double**& fittedPeakInfo, int& fittedPeakCount,
		TApplication* rootApp);
void MyFineFit(int* data, int dataLength, int roughPosition, int roughHeight,
		double**& fittedPeakInfo, int& fittedPeakCount,
		PeakReference* interpolator, int referenceOffset);

int peakCount = 0;

void TCanvas::Closed() {
	Emit("Closed()");

	gApplication->Terminate(0);
}

int main(int argc, char** argv) {
	TApplication* app = new TApplication("app", 0, NULL);

	const int sourceRefSize = noiseRefSize;
	double* sourceRef = new double[sourceRefSize];
	for (int i = 0; i < sourceRefSize; i++)
		sourceRef[i] = noiseReference[i];

	const int numSamples = 100;
	const int numPeaks = 1;
	const int seed = 974488;
	InlineStream IS(numSamples, numPeaks, seed, sourceRef, sourceRefSize);
	int* dataStream = IS.stream;
	double** exactPeakInfo = IS.peakInfo;

	const int peakReferenceSize = referenceSize;//roughGaussSize;
	double* peakReference = new double[peakReferenceSize];
	double peakReferenceMax = 0;
	double peakReferenceSum = 0;
	for (int i = 0; i < peakReferenceSize; i++) {
		peakReference[i] = reference[i];//roughGauss[i];
		peakReferenceSum += peakReference[i];
		if (peakReferenceMax < peakReference[i]) {
			peakReferenceMax = peakReference[i];
		}
	}
	PeakReference* interp = new PeakReference(peakReference, peakReferenceSize);

	double* chisquareArray = new double[numSamples-peakReferenceSize];
	int* positionArrayInt = new int[numSamples-peakReferenceSize];
	double* positionArrayDouble = new double[numSamples-peakReferenceSize];

	//	int histSampleSize = 15; // make an odd number
	//	int deviation = histSampleSize / 2;
	//	TH1I* hist = new TH1I("hist", "hist", histSampleSize, deviation - 0.5,
	//			deviation - .5);

	int fittedPeakCount = 0;
	double** fittedPeakInfo = new double*[fittedPeakCount];

	double chisquare = 10000000000000000;
	double oldChisquare = 10000000000000000;
	double veryOldChisquare = 10000000000000000;
	double height = 0;
	double oldHeight = 0;
	double veryOldHeight = 0;
	int time = 0;
	int oldTime = 0;
	int veryOldTime = 0;
	clock_t startTime = clock();
	for (int i = 0; i < numSamples - peakReferenceSize; i++) {
		veryOldChisquare = oldChisquare;
		veryOldHeight = oldHeight;
		veryOldTime = oldTime;
		oldChisquare = chisquare;
		oldHeight = height;
		oldTime = time;
		chisquare = 10000000000000000;
		double heightGuess;
		double chisquareTmp = 0;

		double sum = 0;
		for (int j = 0; j < peakReferenceSize; j++)
			sum += dataStream[i + j];
		heightGuess = sum / peakReferenceSum;

		double dataSum = 0;
		for (height = heightGuess-heightGuess/10; height < heightGuess+heightGuess/10+1; height += .5) {
			for (int j = 0; j < peakReferenceSize; j++) {
				chisquareTmp += (((height * peakReference[j]) - dataStream[i
						+ j]) * ((height * peakReference[j])
						- dataStream[i + j])) / 7;//(dataStream[i + j]
						//* dataStream[i + j]);
				dataSum+=dataStream[i+j];
			}

			if (chisquareTmp < chisquare) {
				chisquare = chisquareTmp;
			}
		}
		if(dataSum==0){
			dataSum = 0.000000000001;
		}
		chisquare/=fabs(dataSum/peakReferenceSize);

		chisquareArray[i] = chisquare/(peakReferenceSize-3);
		positionArrayInt[i]=i;
		positionArrayDouble[i] = i+8;
		time = i - 6;

//		cout << i + 8 << "\t" << chisquare / (peakReferenceSize - 2) << endl;
		if (oldChisquare < veryOldChisquare && oldChisquare < chisquare) {
			if (chisquare / (peakReferenceSize) < 0.05 && height > 10) {
				if (fittedPeakCount != 0 && !((i - 6
						- fittedPeakInfo[fittedPeakCount - 1][2] >= 14 && i - 6
						- fittedPeakInfo[fittedPeakCount - 1][2] <= 16)
						&& (fittedPeakInfo[fittedPeakCount - 1][0] / height
								>= 24 && fittedPeakInfo[fittedPeakCount - 1][0]
								/ height <= 26))) {
					//if(veryOldChisquare<oldChisquare && newChisquare<oldChisquare){
					//for (int j = -deviation; j < deviation + 1; j++)
					//hist->Fill(j, dataStream[i - j + 4]);
					//RootFineFit(hist, histSampleSize, i+4, height, fittedPeakInfo,
					//	fittedPeakCount, app);
					//hist->Reset();
					addToList(fittedPeakInfo, fittedPeakCount, i - 6, height,
							chisquare);
					//int* fitRegion = new int[peakReferenceSize];
					//for (int j = 0; j < peakReferenceSize; j++)
					//	fitRegion[j] = dataStream[i + j];
					//MyFineFit(fitRegion, peakReferenceSize, i, height,
					//		fittedPeakInfo, fittedPeakCount, interp, 4);
					//}
				} else if (fittedPeakCount == 0) {
					//				if(veryOldChisquare<oldChisquare && goodChisquare<oldChisquare){
					//for (int j = -deviation; j < deviation + 1; j++)
					//hist->Fill(j, dataStream[i - j + 4]);
					//RootFineFit(hist, histSampleSize, i+4, height, fittedPeakInfo,
					//	fittedPeakCount, app);
					//hist->Reset();
					addToList(fittedPeakInfo, fittedPeakCount, i - 6, height,
							chisquare);
					//int* fitRegion = new int[peakReferenceSize];
					//for (int j = 0; j < peakReferenceSize; j++)
					//	fitRegion[j] = dataStream[i + j];
					//MyFineFit(fitRegion, peakReferenceSize, i, height,
					//		fittedPeakInfo, fittedPeakCount, interp, 4);
					//				}
				}
			}
		}
	}
	double elapsedTime = (clock() - startTime) / ((double) CLOCKS_PER_SEC);
	EvaluateAccuracy(exactPeakInfo, numPeaks, fittedPeakInfo, fittedPeakCount);
	//	cout << "elapsed Time = " << elapsedTime << endl;

	TFile f("./output.root","update");
	TGraph g(numSamples-peakReferenceSize,positionArrayInt,dataStream);
	g.Write("stream");
	TGraph g2(numSamples-peakReferenceSize,positionArrayDouble,chisquareArray);
	g2.Write("chisquare");
	ofstream fitfile("fittedPeaks");
	for (int i = 0; i < fittedPeakCount; i++)
		fitfile << fittedPeakInfo[i][2] << "\t" << fittedPeakInfo[i][0] << "\t"
				<< fittedPeakInfo[i][4] / (peakReferenceSize - 2) << "\t"
				<< fittedPeakInfo[i][5] << endl;
	fitfile.close();

	ofstream peakfile("peaks");
	for (int i = 0; i < numPeaks; i++)
		peakfile << exactPeakInfo[i][1] << "\t" << exactPeakInfo[i][0] << endl;
	peakfile.close();

	ofstream streamfile("stream");
	for (int i = 0; i < numSamples; i++)
		streamfile << i << "\t" << dataStream[i] << endl;
	streamfile.close();

	delete interp;

	return 0;
}

void addToList(double**& list, int& count, double time, double height,
		double chi2) {
	peakCount++;
		if (peakCount == 1)
			cout << peakCount << " peak identified\r" << flush;
		else
			cout << peakCount << " peaks identified\r" << flush;

	double** newList = new double*[count + 1];
	for (int i = 0; i < count; i++)
		newList[i] = list[i];
	delete list;
	newList[count] = new double[6];
	newList[count][0] = height;
	//	newList[count][1] = heightError;
	newList[count][2] = time;
	//	newList[count][3] = timeError;
	newList[count][4] = chi2;
	newList[count][5] = 0;
	count++;
	list = newList;
}

int numPeaksRemoved;

void RemoveFromList(double**& list, int& count, int indexToRemove) {
	double** newList = new double*[count - 1];
	for (int i = 0; i < indexToRemove; i++)
		newList[i] = list[i];
	for (int i = indexToRemove; i < count - 1; i++)
		newList[i] = list[i + 1];
	count--;
	numPeaksRemoved++;
	list = newList;
}

void MyFineFit(int* data, int dataLength, int roughPosition, int roughHeight,
		double**& fittedPeakInfo, int& fittedPeakCount, PeakReference* interp,
		int referenceOffset) {

	double chisquare = 10000;
	double chisquareTmp = 0;
	double finePosition = 0;
	int height = 0;
	for (int i = roughHeight - 10; i < roughHeight + 11; i++) {
		for (double j = -1; j <= 1; j += 0.01) {
			chisquareTmp = 0;

			for (int k = 0; k < dataLength; k++)
				chisquareTmp += ((i * interp->cubicInterpolationAt(
						referenceOffset + k - j) - data[k]) * (i
						* interp->cubicInterpolationAt(referenceOffset + k - j)
						- data[k])) / (data[k] * data[k]);

			if (fabs(chisquareTmp - 1) < fabs(chisquare - 1)) {
				chisquare = chisquareTmp;
				finePosition = j;
				height = i;
			}
		}
	}

	if (chisquare / (dataLength - 3) < 1.01)
		addToList(fittedPeakInfo, fittedPeakCount,
				roughPosition + finePosition, height, chisquare);
}

void RootFineFit(TH1I* hist, int histSampleSize, int roughPosition,
		double roughHeight, double**& fittedPeakInfo, int& fittedPeakCount,
		TApplication* app) {

	TF1* myGaus = new TF1("myGaus", "gaus(0)", -histSampleSize + 100,
			histSampleSize - 100);
	myGaus->FixParameter(2, 2.3887175);
	myGaus->SetParLimits(2, 0, 100);
	myGaus->SetParameter(2, 2.3887175);
	myGaus->SetParameter(0, roughHeight);
	myGaus->SetParLimits(0, 0, 16384);
	myGaus->SetParLimits(1, -1, 1);
	myGaus->SetParameter(1, 0);

	TFitResultPtr fit = hist->Fit("myGaus", "MSQB0");

	double chi2 = fit->Chi2();
	//	cout << chi2 << endl;
	double height = fit->Value(0);
	double offset = fit->Value(1);

	if (chi2 / (histSampleSize - 4) < height / 100) {
		//		app->Run(true);
		addToList(fittedPeakInfo, fittedPeakCount, (double) roughPosition
				+ fit->Value(1), fit->Value(0), chi2);
	}
}

int numPeaksMissed = 0;
int numFalsePeaks = 0;

void EvaluateAccuracy(double** exactPeakInfo, int exactPeakCount,
		double**& fittedPeakInfo, int& fittedPeakCount) {
	cout << endl;
	double** accuracy = new double*[exactPeakCount];
	for (int i = 0; i < exactPeakCount; i++)
		accuracy[i] = NULL;

	for (int i = 0; i < fittedPeakCount; i++) {
		double sigma = 100000000000000;
		double sigmaTemp = 0;
		int index = 0;

		for (int j = 0; j < exactPeakCount; j++) {
			sigmaTemp = fittedPeakInfo[i][2] - exactPeakInfo[j][1];

			if (fabs(sigmaTemp) < fabs(sigma && fabs(sigmaTemp) < 1)) {
				sigma = sigmaTemp;
				index = j;
			}
		}

		if (accuracy[index] == NULL) {
			accuracy[index] = new double[3];
			accuracy[index][0] = sigma;
			accuracy[index][1] = (fittedPeakInfo[i][0]
					- exactPeakInfo[index][0]) / exactPeakInfo[index][0];
			accuracy[index][2] = i;
			fittedPeakInfo[i][5] = 0;
		} else if (fabs(sigma) < fabs(accuracy[index][0])) {
			fittedPeakInfo[(int) accuracy[index][2]][5] = 1;
			accuracy[index][0] = sigma;
			accuracy[index][1] = (fittedPeakInfo[i][0]
					- exactPeakInfo[index][0]) / exactPeakInfo[index][0];
			accuracy[index][2] = i;
			fittedPeakInfo[i][5] = 0;
		} else {
			fittedPeakInfo[i][5] = 1;
		}
	}

	TFile f("./accuracy.root","recreate");
	TH1I histo("hist","hist",128,-0.5,16348-0.5);
	ofstream accuracyfile("accuracy");
	for (int i = 0; i < exactPeakCount; i++) {
		if (accuracy[i] != NULL) {
			accuracyfile << i << "\t" << accuracy[i][0] << "\t"
					<< accuracy[i][1] << "\t" << accuracy[i][2] << endl;
		} else{
			histo.Fill((int)exactPeakInfo[i][0]);
			numPeaksMissed++;
		}
	}
	histo.Write("accuracyHist");
	accuracyfile.close();

	for (int i = 0; i < fittedPeakCount; i++) {
		if (fittedPeakInfo[i][5] != 0)
			numFalsePeaks++;
	}
	cout << numFalsePeaks << " false peaks\t" << numPeaksMissed
			<< " peaks missed" << endl;
	/*
	 if (numPeaksRemoved != 1)
	 cout << numPeaksRemoved << " peaks removed" << endl;
	 else
	 cout << numPeaksRemoved << " peak removed" << endl;
	 if (fittedPeakCount != 1)
	 cout << fittedPeakCount << " peaks remaining" << endl;
	 else
	 cout << fittedPeakCount << " peak remaining" << endl;
	 if (numPeaksMissed != 1)
	 cout << numPeaksMissed << " peaks missed" << endl;
	 else
	 cout << numPeaksMissed << " peak missed" << endl;
	 if (numFalsePeaks != 1)
	 cout << numFalsePeaks << " false peaks" << endl;
	 else
	 cout << numFalsePeaks << " false peak" << endl;
	 */
}
